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Abstract 

We present an approximate deconvolution (AD) large eddy simulation (LES) model for the two-layer quasi- 
geostrophic equations. We applied the AD-LES model to mid-latitude two-layer square oceanic basins, which are 
standard prototypes of more realistic stratified ocean dynamics models. Two spatial filters were investigated in the 
AD-LES model: a tridiagonal filter and an elliptic diff'erential filter. A sensitivity analysis of the AD-LES results with 
respect to changes in modeling parameters was performed. The results demonstrate that the AD-LES model used in 
conjunction with the tridiagonal or diff'erential filters provides additional dissipation to the system, allowing the use 
of a smaller eddy viscosity coefficient. Changing the spatial filter makes a significant difference in characterizing the 
effective dissipation in the model. It was found that the tridiagonal filter introduces the least amount of numerical 
dissipation into the AD-LES model. The differential filter, however, added a significant amount of numerical dissipa- 
tion to the AD-LES model for large values of the filter width. All AD-LES models reproduced the DNS results at a 
fraction of the cost within a reasonable level of accuracy. 

Keywords: Approximate deconvolution; Large eddy simulation; Subfilter- scale parameterization; Two-layer 
quasigeostrophic equations; Forced-dissipative ocean models; Large-scale ocean circulation. 



1. Introduction 



The investigation of characteristics of forced-dissipative general circulation models is of primary importance in 
developing our understanding of the large-scale nonlinear motions of geophysical flows. As one of the main circulation 
sources, winds drive the general circulation associated with the subtropical and subpolar gyres, which can be identified 
with the strong, persistent, sub-tropical and sub-polar western boundary currents in the North Atlantic Ocean (the Gulf 
Stream and the Labrador Current) and North Pacific Ocean (the Kuroshio and the Oyashio Currents) and sub-tropical 
counterparts in the southern hemisphere (Stommel 1972[ Mc Williams 2Q06| ). One of the major similarities between 
the various ocean basins is the asymmetry of the gyres: strong western boundary currents and weaker flow in the 
interior; weak and shallow eastern boundary currents. The most obvious motivation for being interested in forced- 



dissipative wind-driven ocean circulation is the connection between ocean currents and climate dynamics ( |Ghil et aL 
[2008] ). 

The wind-driven circulation in an enclosed, midlatitude rectangular or square basin is a classical problem, studied 
extensively by modelers ( |AlIen| p^SOt |HoUand and Rhines| p^80l [Griffa and Salmon| p^89l [Wllsl [20061 IMifler 



2007| ). Various models are derived from the full-fledged equations of geophysical flows, Boussinesq equations (BEs) 



or the primitive equations (PEs), to guide the theoretical studies on boundary currents, alternating zonal flows, or jet 
formations, as well as to identify some key issues related to the relative insensitivity of the model dynamics to the 
changes of parameters that is closely linked to a dynamical system point of view ( [Speich et"aL 1 995 [ |Meacham| [2000 



[Chang et al.) [2QQT| [Nauw et al.[ [20041 [Dijkstra| [2005] [Dijkstra and Ghil|[2QQ5] ). The quasigeostropic (QG) model is 
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a simplification of the primitive equation model that retains many of the essential features of geophysical fluid flows. 
Details of the mathematical and physical approximations may be found in standard textbooks on geophysical fluid 



dynamics, such as Pedlo sky, ( |1987| ), |Vallis| ( |2006 ), and McWilliains| ( [2006) ). The main assumptions that go into the 
QG models are: the hydrostatic balance, the yS-plane approximation, the geostrophic balance, and the eddy viscosity 
parameterization . 

The one-layer QG model, sometimes called the barotropic vorticity equation (BVE), represents one of the most 
commonly used mathematical models for these types of geostrophic flows with various dissipative and forcing terms 
( Majda and Wang| 2QQ6t Vallis 2006[ Nadiga and Margolml 2001 ). In reality, the ocean is a stratified fluid on a 
rotating Earth driven from its upper surface by patterns of momentum and buoyancy fluxes ( [Marshall et al.[pW7] ). 
While the barotropic model is not stratified, it exhibits many of the features that are observed in the stratified case. 
To explore some of the efl'ects of the stratification, the one-lay er barotropic equation can be extended to the 1.5-layer 
model, also called the reduced gravity QG model ( Qzgokmen et al. 200l| . There are two layers in this model, but 
the second layer is infinitely deep and at rest (passive), and the dynamics are eff'ectively barotropic. The two-layer 
model takes the next step in incre asing the complexity of stratification by addin g a seco nd dyna mically active layer 
Holland||l978[|Qzgokmen and Chassignet|| 1998, Berloff^ and McWiniams,.1999|]DiBattista and M^da| [2001 [[Berloff 



et al. 2009| ). The dynamics in this model include the first baroclinic modes. The complexity of the models could be 
increased by adding more active layers, resulting in the N-layer models (Siegel et aLl[2001| ), which, in turn, yield the 
three dimensional primitive equations when N goes to infinity (Mc Williams,, 2006) ). In this study, we use the two-layer 
QG (QG2) model. 

Geophysical turbulence is strongly aff"ected by the planetary vorticity, the variation of the Coriolis parameter with 
latitude, the so-called effect ( [Maltrud and ValHsl \TWl\ [Smith et aT] [20021 [Chen et al.[ [2003] ). The inverse cascade 
typically occurring in pure two-dimensional turbulence, in this case preferentially transfers small-scale energy towards 
zonal modes; the resulting flow is then anisotropic and characterized by a strong interaction between waves and 
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42 turbulence, and is known as the arrest of the inverse energy cascade ( jRhines[ [1975[ [Sukoriansky et al.[ [2007 [ [Espa 



et al.[[2008{[San and Staples[[2013b[ ). [Rhines[ ( [l975[ ) explained the emergence of flow anisotropy and the organization 



of a banded pattern of alternating zonal currents, or jets, due to Rossby wave dynamics in terms of a competition 
between nonlinear and 13 terms in the barotropic vorticity equation. Under the effects of planetary rotation, Rossby 
waves dominate turbulent motions prohibiting the triad interactions, and arrest the inverse energy cascade when the 
scale of motions becomes larger than a critical value, later known as the Rhines scale (Tanaka an d Akitomo '2010 ). 

Along with the Rhines scale which is a measure of the strength of nonlinear interactions, another important 
scale for determining the dynamics of the large scale motions in the ocean is the Munk scale ( [Munk[ [1950[ ), which 
corresponds to the dissipative behavior of the system and can be linked to the Reynolds number. Although the water 
molecular viscosity is around 10"^ m^s~^, the one- and two-layer QG models use viscosities on the order of 10^ m^^"^ 
This is called eddy viscosity (EV) parameterization, and is used because the horizontal scale of the ocean basin is 
much larger than the effective scale for molecular diffusion. An impractically fine resolution would be necessary if 
the ocean models were to resolve the full spectra of turbulence down to the Kolmogorov scale. Thus, the viscosity 
coefl&cients employed in the QG models typically remain much greater than the molecular viscosity ( [Campin et al. 



2011 ). The eddy viscosities generally used in the oceanic models are summarized in Table [T] The eddy viscosity 



iparameteriza tion used in the QG models plays a crucial role in the dynamics of the problem. Indeed, [Berloff and 
[McWilliams (1999]) studied the wind-driven circulation in a three-layer QG model for varying values of the eddy 
viscosity coefficient in a square oceanic basin. For y = 1200 m^s~^ an asymmetric steady state was found. When the 
eddy viscosity coefficient was decreased, the flow first displayed a variability characterized by the presence of interior 
Rossby waves. At v = 1000 m^s~^, the flow regime showed a quasi-periodic variability. At a smaller eddy viscosity 
coefficient, starting from y = 800 m^^"^ the flow regime was chaotic and showed a persistent eastward jet penetration 
by fluctuating between two preferred states, one of which corresponds to a low energy state and a long eastward jet, 
and the other to a high energy state and a short jet. The study of Berloff and McWilliams[ ( [l999[ ) clearly shows that 
different EV coefficients can result in different dynamics of the QG models. Thus, a natural question is "What EV 
coefficient should be used in the QG models?" The EV coefficients summarized in Table [T] seem to convey, at first 
glance, a confusing message: they vary by as much as an order of magnitude. At a closer look, however. Table [T] 
clarifies this issue: With the ever increasing computational power, the mesh size used in numerical simulations with 
the QG models constantly decreases and allows the use of smaller EV coefficients. The development of a rigorous, 
mathematical understanding and subsequent modeling strategy for the eddy viscosity coefficients (see Table [T]) is the 



'elephant in the room," one of the major unsolved problems in ocean model ing ([Visbeck et al.[p^97{|Campin et al. 
2011[ |Majda and Wang[ |2006[ [Cushman-Roisin and Beckers| |2009t |Valhs| |2006| ). Although addressing this grand 



challenge is beyond the scope of this report, we do address the intimate relationship between the EV coefficients and 
the numerical resolution employed by the QG models. 



Study 



Table 1 : The eddy viscosity coefficients used in QG models 

Range of y (m^^"^) Resolution 



Hollar 




Jiang et al. '{ 1995 



(1975) 



Qzgokmen an d Chassign et 



Berloff and Mc\ 
Sura etal.,(,200T| 



391 



Tanaka and Akitomo (20101 



10000 
- 10000 



500- 
6000 
330 
300 
50 

400 - 1600 
200 
100 
100 



40x80 

74x50 

50x50 

50x100 

151x151 

256x256 

120x120 

512x256 

500x500 
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79 
80 
81 
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102 
103 
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To capture the under-resolved flow, i.e., the flow in the regions where the grid size becomes greater than the 
specified Munk scale, large eddy simulation (LES) appears as a natural choice. Most of the LES models have been 



developed for three-dimensional turbulent flows, such as those encountered in engineering applications (|Sagaut[[2006 



Berselli et al.[|2006] ). These LES models fundamentally rely on the concept of the forward energy cascade and so their 



extension to geophysical flows is beset with difficulties. The efl'ective viscosity values in oceanic models are much 
greater than the molecular viscosity of seawater, hence a uniform eddy viscosity coefficient is generally used to 
parameterize the unresolved, subfilter- scale eff'ects in most oceanic models ( Mc Williams | 2006 Vallis 2006| ). LES 
models specifica lly developed for two-dimensional turbulent flows, such as those in the ocean and atmosphere, ar e 



83 relatively scarce ( Fox-Kemper and Menemenlis 



2008 



Awad et al. 



2009 



Ozgokmen et al. 



2009 



Chen et al. 



at least when compared to the plethora of LES models developed for three-dimensional turbulent flows. Ho 



20111, 



m and] 



|Nadiga| ( |2003| ) combined the uniform eddy viscosity parameterization with the alpha regularization LES approach to 
capture the under-resolved flow where the grid length becomes greater than the specified Munk scale of the problem. 
In that work, the structural alpha parameterization was tested on the barotropic vorticity equation (BVE) in an ocean 
basin with double-gyre wind forcing, which displays a four-gyre mean ocean circulation pattern. It was found that the 
alpha models provide a promising approach to LES closure modeling of the barotropic ocean circulation by predicting 
the correct four-gyre circulation structure for under-resolved flows. 

San et al.|p011| ) put forth a new LES closure modeling strategy for two-dimensional turbulent geophysical flows. 



The new closure modeling approach utilizes approximate deconvolution (AD), which is particularly appealing for 
geophysical flows because of no additional phenomenological approximations to the BVE. Similar to the method 
suggested by |Holm and Nadiga| ( |2003| ), this framework also uses a Laplacian operator with a constant eddy viscosity 
coefficient to account for the dissipation mechanism. For a given system with eddy viscosity dissipation, the subfilter- 
scale contribution, however, is modeled by a non eddy viscosity AD closure approach. The AD approach can achieve 
high accuracy by employing repeated filtering, which is computationally efficient and easy to implement. The AD 



98 method has been used successfully in LES of three-dimensional turbulent engineering flows ( Stolz and Adams | 1999 



|Stolz et aL| |2001a|b[ |2004[ [Domaradzki and Adams 2002 ) and even of small scale geophysical flows, such as the 
atmospheric boundary layerl Chow et al. , 2005; Chow and Street| [20091 |Duan al.|[20TQl[Zhou and Chow|[2QTT] ). 
The AD methodology was also used in LES of large scale geophysical flows, such as the barotropic ocean circulation 
flow. To assess the new AD closure modeling approach, San et al." ("2011) tested it on the same two-dimensional 
barotropic flow problem as that employed in |Nadiga and Margolin ( 2001) and in|Holm and Nadig"a|p003). It was 



shown that the new LES -AD model provides an accurate approximation for under-resolved subfilter- scale effects. 



The main goal of this report is to extend the LES-AD approach used for the one-layer QG model ( San et al. 201 1 



to the two-layer QG model. A quantitative analysis of the eff'ects of using the AD-LES model on QG2 models was 
performed in conjunction with the tridiagonal and diff'erential filters. We investigated whether the combination of 
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LES-AD modeling and a particular spatial filter can, in fact, account for some of the eddy viscosity parameterization 
used in practical QG numerical simulations. Our numerical experiments show that the AD-LES model does add 
numerical dissipation, but the exact amount and form still need to be determined. A sensitivity analysis was performed 
to find out how much of the dissipation the AD-LES model, equipped with various spatial filters, can account for. We 
demonstrated that the amount of the dissipation added to the system depends on the free modeling parameters. We 
emphasize that this issue is common to LES modeling in general. Indeed, not only is it hard to find the "best" LES 
model, i.e., the model that produces the most accurate results at the lowest computational cost, but once this model is 
found, it is often hard to decide whether the success of the model is due to the actual closure model or the numerical 
discretization used ( Berselli et al.[|20Q6{|Sagaut[|2006) . In an actual LES of turbulent flow there are several ingredients 



some are used at the continuum level (e.g., the closure model with its various parameters), and some are used at 
the discrete level (e.g., the temporal and spatial discretization or the linear solver). Often, it is hard to disentangle the 
modeling efl'ects from the numerical discretization efl'ects. Our QG setting is no difl'erent in this respect. We plan to 
investigate this complex relationship in a future study, by performing extensive numerical experiments in simplified 
settings and by developing mathematical support for these numerical results. 

The rest of the paper is organized as follows: Section |2] presents the two-layer QG equations for large-scale 
geophysical flows. The proposed AD methodology, which yields the mathematical model used in this report, is 
presented in Section [3] The numerical methods used in our simulations are briefly discussed in SectionH] The results 

125 for the new AD model are presented in Section[5] Finally, the conclusions are summarized in Section|6] 

126 2. Governing equations 

127 2.7. The two-layer quasigeostrophic equations 

128 The two-layer quasigeostrophic model used in this study is one of the simplified forced-dissipative oceanic models 

129 that considers baroclinic eff'ects. The stratified ocean is partitioned into two isopycnal layers, each of constant depth, 

130 density and temperature. The governing quasigeostrophic potential vorticity equations for the two dynamically active 

131 layers are ( |Pedlosky| pW/j |Salmon| p^98[ |Mc Williams [ 12506) ) 



^+/(<Ai,^i) = Di+Fu (1) 
ot 

^-^+J{i^i.qi) = + (2) 
ot 

where the layer index starts from top, qt represents potential vorticities, and ij/i denotes for streamf unctions. The Jaco- 
bian operator is defined as /(a, ^) = |f ^ - ^ • The dissipation and forcing (Ekman pumping) terms are represented 
by Di, and F/, respectively. The potential vorticities for each layer are related to the velocity streamfunctions through 
the following elliptic coupled system of equations: 

qi = vVi+ySj+-7r(^2-^i), (3) 



f 

qi = VV2+ySj+^(<Ai-^2). (4) 

136 The isopycnal flow velocity components can be found from the velocity streamfunctions: 

oy ox 

137 The two symbols /3 and /o are parts of the linearized yS-plane approximation to the Coriolis parameter / = /q -h j3y. 

138 Here /o = 2 D sin(0o) is the local rotation rate aty = 0, where D is the rotational speed of the earth and 0o is the 

139 latitude at 3; = 0. This is equivalent to approximating the spherical Earth with a tangent plane at 3; = 0. Stratification 

140 is represented by two stacked isopycnal layers with thicknesses Hi and H2, starting from the top, and g' = 

141 is reduced gravity associated with the density jump between the two layers in which Ap is the density diff'erence 

142 between the two layers, pi is the reference (upper layer) density, and g is the gravitational acceleration. The inertial 
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143 radius of deformation between layers, a measure of stratification strength, is defined as the Rossby deformation radius 

144 Rd = -yj ^'fi^^ ^ where H = H\ + H2. In this study, the top and bottom layers of the ocean are forced by an Ekman 

145 pumping of the form 

Fi = -^k-Vxf, (6) 

F2 = -rvV2, (7) 

146 where r = (t^^\ t^^^) is the stress vector for surface wind forcing, and k is unit vector in vertical direction. In the present 

147 model, we use a double-gyre wind forcing only for zonal direction: r^^^ = tq cos^ ^yj^ where L is the meridional length 

148 of the ocean basin centered at j = 0, and tq is the maximum amplitude of the wind stress. This form of wind stress 

149 represents the meridional profile of easterly trade winds, mid-latitude westerlies, and polar easterlies from South to 

150 North. The bottom Ekman layer is parameterized by a linear bottom friction with coefficient y. In the equations 

151 above, V and are the gradient and Laplacian operators, respectively. For the dissipation terms, the following EV 

152 parameterizations are used: 

D, = vV^i/iu (8) 
D2 = VVV2, (9) 

153 where v is eddy viscosity coefficient. 

154 2.2. Governing equations in dimensionless form 

155 The governing equations can be written in dimensionless form by using the Sverdrup balance to set the velocity 

156 scale of the form 

V=^^. (10) 

157 The dimensionless variables (denoted by tilde) are defined as 

158 Then the two-layer quasigeostrophic equations in dimensionless form become 

^^J(iffuqi) = D,^sm(2nyX (12) 
at 

^^Jih.h) = D2-crV^2, (13) 
ot 



159 in which the dissipative terms can be written as 



Di = AVVi (14) 



D2 = AV4^2. (15) 

160 In dimensionless form, the kinematic relationships between potential vorticities and streamfunctions become: 

Fr 

qi = RoV^i/^i +5;+ — ((A2-1A1), (16) 
o 

q2 = RoV2(/r2+^+ 7^(iAl -<A2). (17) 

1-0 

161 For clarity of exposition, in the remainder of the paper we will drop the tilde symbol used for the dimensionless 

162 variables. In the two-layer QG model, S = ^ is the aspect ratio of vertical layer thicknesses, Ro is the Rossby 
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163 number, Fr is the Froude number, A is the lateral eddy viscosity coefficient, and cr is the Ekman bottom later friction 

164 coefficient. The definitions of these dimensionless parameters are: 
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V fnV VL y y 

Ro = — ; Fr = Re = — ; A = — ; o- = (18) 

The following three length scales are useful for setting the problem parameters: (i) the Munk scale, = ( |j , 
for the viscous boundary layer; this is related to the smaller scale dissipation; (ii) the Stommel scale, Ss = ^, for the 

bottom boundary layer thickness; this is accounting for larger scale damping; and (iii) the Rhines scale, ^/ = I ^ J , 

for the inertial boundary layer; this is measuring the strength of the nonlinearity. 

In order to complete the mathematical model, boundary and initial conditions should be prescribed. In many 
theoretical studies of ocean circu l ation, t he modelers either use free- slip bo undary conditions or no- slip boundary 
conditions. Following |Cummins| (|l992|; [Oz gokmen and Chassign et] (|l998|, we use free- slip boundary conditions 
for the velocity for both isopycnal layers, which translates into homogenous Dirichlet boundary conditions for the 
vorticity (Laplacian of streamf unction): V^</^|q = 0. The impermeability boundary condition is imposed as = 0. 
We start from a rest state (i// = 0), integrate the model until a statistically steady state is obtained, and continue for 
several decades to compute time-averaged results. 



176 3. Approximate deconvolution method 

177 The goal in AD is to use repeated filtering in order to obtain approximations of the unfiltered unresolved flow 

178 variables when approximations of the filtered resolved flow variables are available. These approximations of the 

179 unfiltered flow variables are then used in the subfilter- scale terms to close the LES system. To derive the new AD 

180 model, we start by denoting by G the spatial filtering operator: Gu = u, Gu = u and so on, where u represents any 

181 flow variable (i.e., potential vorticity and the streamfunction in this study). Since G = I - (I - G), an inverse to G can 

182 be written formally as the non-convergent Neumann series: 

oo 

G-'-Yj^I-Gy. (19) 

183 Truncating the series gives the van Cittert approximate deconvolution operator, Q^. We truncate the series at N and 

184 obtain as an approximation of G"^ : 

QN = Yj^i-Gr\ (20) 

185 where / is the identity operator. The approximations are not convergent as N goes to infinity, but rather are 



186 asymptotic as the filter radius. A, approaches zero ( [Berselli et al.| [2006 J . An approximate deconvolution of any 

187 variable u can now be obtained as follows: 

t/* = Qnu. (21) 

188 For higher values of N, we get increasingly more accurate approximations of u\ 

Qi = / (22) 

Qi = 21 -G (23) 

23 = 3/-3G-fG^ (24) 

24 = 4/-6G-h4G^ -G^ (25) 

25 = 5/ - lOG -h lOG^ - 5G^ -h G"^ (26) 
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Following the same approach as that used in |Dunca and Epshteyn| ( |2006] ), one can prove that these models are 



190 highly accurate (0(A^^^^) modeling consistency error) and stable. For example, if we choose N = 5,we can find an 

191 AD approximation of the resolved variable q as 

q* = 5q- lOq + 10§ - 5| + | (27) 

192 and, similarly, an AD approximation of the variable i// as 

i//* = 5i^- 10</^ + 10<^ -5^ + ^, (28) 

193 where q and i// are the resolved potential vorticity and streamf unction variables. We use a bar to denote the application 

194 of one filtering operation. Using ( [27] ) and ( [28] ), we can now approximate the subfilter- scale contribution by applying 

195 a filter to the governing equation. This results in the following model: 

% + /(iAi,^i) = AVVi + sin(27r};) + ^*, (29) 
ot 

§ + /(fe^2) = AVV2-c^vV2+^;, (30) 
ot 

196 where S* is the subfilter- scale term for the layer, given by 



^; = -/((A;,<) + /(^/,^/), (31) 

197 where asterisk represents the approximated value for the unfiltered (unresolved) quantities. To completely specify the 

198 new AD model ([29|)-([3T|), we need to choose a computationally efficient filtering operator. In Section |5] we will show 

199 that the selection of the filtering operator aff'ects the dissipative behavior of the system. 

200 3.1. Tridiagonal filter 

201 Following [Stolz and Adams ( 1999 ), we use the following discrete second-order tridiagonal filter (TF): 



where fi represents the filtered value of a discrete quantity fi. Here, the subscript / is the spatial index in the x- 
direction. This results in a tridiagonal system of equations for each fixed value of y. After solving Eq. ( [32| ), we use the 
same filter in the }^-direction (i.e., we replace index / with j) for each fixed value of x. The resulting tridiagonal system 
of equations is solved efl&ciently by using the well-known Thomas algorithm. Since the TF has been constructed in the 
physical space, a Fourier analysis is applied to study its characteristics in the wavenumber space. This analysis leads 
to the transfer function, T(oj), that correlates the Fourier coefficients of the filtered variable to those of the unfiltered 
variable as follows: 

fk = T(aj)fk, (33) 

where fk and fk are the Fourier coefficients of the filtered and unfiltered variables, respectively (i.e., /i = 2 fk^^^^' and 
fi = Yj fk^^^^'^ where xi = AJ and A;^ is the grid spacing in the x-direction). Using the relation cos(^) = (e^^ -h e~^^)/2, 
the transfer function of the TF given in Eq. ([32]) can be written as 



r<->(.) = (i..)-lt^, (34) 

\ 2 I 1 -\-2a cos((jL)) 

where oj = kAx is the modified wavenumber in the x-direction. The free parameter, a, which is in the range < < 
0.5, determines the filtering properties, with high values of a yielding less dissipative results. If the transfer function 
of the filter used in the AD closure is positive, then the existence and uniqueness of strong solutions of the AD model 
can be proved ( Stanculescu| 2008] ). The transfer function corresponding to the TF becomes positive definite in the 



interval of < \a\ < 0.5. More details can be found in |San et aL] ( |2011| ). 



7 



1.2 I- 



0.8 - 



0.6 - 



3^ 



0.4 - 



0.2 - 



- 



X, 



^ 



Fourier cut-off filter 



■ TF:a = 0.45 

TF:a = 0.35 

TF:a = 0.25 

TF:a = 0.15 

TF:a = 0.05 



\ 



\ 



\\\ ^ 



\ \ \ \ » 



-0.2, 



J_ 



0.2 0.4 0.6 

co/tc 



0.8 



Figure 1: Transfer functions of the TF for different values of the parameter a. The transfer function of the Fourier cut-off filter is also included for 
comparison purposes. 

To show the characteristics of the TF in Eq. ([32]), we plot in Fig.[l]its transfer function T^^^\(jS) (which is given by 
Eq. ([34])) for different values of the free parameter a. The transfer function of the Fourier cut-off filter is also shown 



219 for comparison purposes (see |Najjarand Tafti (1996^ ). It is known that the Fourier cut-off filter removes the small 

220 scales with wavenumbers cdIti > 1/2, while retaining the larger scales with wavenumbers coin < 1/2. It is clear from 

221 Fig. [T] and Eq. ( [34] ) that a plays the role of a cut-off wavenumber for the TF: a = 0.5 turns off the filter, whereas low 

222 a values result in more dissipation (i.e., high attenuation of all the wavenumber components). 

223 3.2. Elliptic differential filter 

224 The second filter used in our numerical investigation is the elliptic differential filter (DF) ( Germano 1986[ Sagaut 

225 [20061 IBerselli et al.|[2QQ6] ): 



f = f on5Q, (36) 

where Q is the computational domain and A is the Helmholtz length, which determines the effective width of the filter. 
The DF is also called Helmholtz filter. The filtered value / is obtained by applying the inverse Helmholtz operator to 
the unfiltered flow variable /. This inversion is done efficiently by using the fast Fourier transform (EFT) techniques 
( Press et aLl|1992| ). Specifically, we use the fast sine transform to solve the discrete version of Eq.([35]), which can be 



written as follows: 



fij - ^ I Kl ^ Af 1 = ^^^^ 
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Figure 2: Transfer functions of the DF for different values of the Helmholtz length A. The transfer function of the Fourier cut-off filter is also 
included for comparison purposes. 



231 The two-dimensional form of the DF in Eq. ( [37] ) is used throughout the paper. In this section, however, to study 

232 the characteristics of the DF in the wavenumber space, we consider the one-dimensional version of the DF (in the 

233 x-direction) 



236 
237 
238 
239 
240 
241 
242 
243 
244 
245 
246 
247 



2 I fi+l ~ '^fi + fi-l 



A? 



(38) 



234 and perform a Fourier analysis similar to the analysis presented in Section 3.1 Thus, the transfer function of the DF 

235 becomes 

T^'''\^) = „ V . ,. • (39) 



1 



|j(2cosM-2) 

It is obvious that the transfer function T^^^^ in Eq. ([39]) is positive, which ensures the well-posedness of the AD 
model (see Stanculescu, (2008 ) ). To show the characteristics of the DF, we plot in Fig. [2] its transfer function, T^^^\ 
for different values of the parameter A. In this study, we parameterize the Helmholtz length /I as a linear function 
of the grid spacing h = Ax = Ay. Thus, increasing the value of A in Fig. |2] amounts to increasing the filter width, 
while keeping the grid spacing fixed. Fig. [2] clearly shows that increasing A (i.e. increasing the filter width) results 
in a significant increase of the dissipation of the DF (the attenuation of the wavenumber components of the filtered 
variable). 

The DF ([35])-([36]) was introduced in LES by |Germano| ([19861). Since then, it has been success fully used in LES 



of three-dimensional engineering flows (Iliescu and Fischer 2003) and small scale oceanic flows ( Qzgokmen et al. 
^009 ). It has also been analyzed mathematically (Dunca and Epshteyn 2006 ; Layton and Lewando wski[ 200 6 ; Layton[ 



and Neda|[2007[[L"ayton and Rebholz||2012[|sranculescU||2008,,Berselli and Lewandowski^|2011j ). In this study, the 
DF is used in the LES of large scale oceanic flows. 
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248 4. Numerical methods 



249 In many physically relevant situations, where the Munk and Rhines scales being close to each other, the solutions 

250 to oceanic models, such as the QG2 models, do not converge to a steady state as time goes to infinity f Medjo|[2000| ). 

251 Rather they remain time dependent by producing statistically steady state with one or multiple equilibria. Therefore, 

252 numerical schemes designed for numerical integration of such phenomena should be suited for such behavior of the 

253 solutions and for the long-time integration. In this study, the governing equations are solved by a fully conservative 

254 finite diff'erence scheme along with a third-order Runge-Kutta adaptive time stepping algorithm. An efficient, linear- 

255 cost, fast sine transform method is utilized for solving the linear coupled inversion subproblem. 

256 4.1. Arakawa scheme for the Jacobian 

257 |Arakawa| (|1966|) suggested that the conservation of energy, enstrophy, and skew- symmetry is sufliicient to avoid 



258 computational instabilities stemming from nonlinear interactions. The second-order Arakawa scheme for the Jacobian 

259 (thenonlinear term in the governing equations) is 

/((A, q) = q) + hii/^, q) + U^lf, q)\ (40) 

260 where the discrete Jacobians have the following forms: 

q) = ^ [-(qi+ij - qi-ijMij+i - ^uj-i) 

Hqij+i - qij-i)(i//i+ij - (41) 
/2(<A. q) = ^^^^ [-qi+i,0i+i,j+i - ^i+ij-i) + qi-i,Mi-i,j+i - ^i-ij-i) 

-\-qij+i(i//i+ij+i - - qij-i(i//i+ij-i - (42) 

Jsil//, q) = -^^-^[-qi+ij+iiif/ij+i - + qi-ij-i(i//i-ij - i/^ij-i) 

-\-qi-ij+i(i//ij+i - i/^i-ij) - qi+ij-i(i//i+ij - i/^ij-i)]- (43) 

261 Note that Ji , which corresponds to the central second-order diff'erence scheme, is not sufl&cient for the conservation 

262 of energy, enstrophy, and skew-symmetry by the numerical discretization. |Arakawa| ( [1966| ) showed that the judicious 

263 combination of /i, /2, and J3 in Eq. (|4Q|) achieves the above discrete conservation properties. 



264 4.2. Time integration scheme 

265 For the time discretization, we employ an optimal third-order total variation diminishing Runge-Kutta (TVDRK3) 

266 scheme ( |Gottlieb and Shu^|1998, ). For clarity of notation, we rewrite the governing equations in the following form: 

t=^" (44) 
at 

267 where subscript / represents the layer index and Ri denotes the discrete spatial derivative operator, including the 

268 nonlinear Jacobian of the convective term, the linear biharmonic diff'usive term, the forcing term, and the subfilter- 

269 scale term. For each layer, the TVDRK3 scheme then becomes: 

= q^^AtRf, 

^(2) = + i^d) + iA,7?^i), (45) 

^' 3 ' 3 • 3 

270 where At is the adaptive time step size, which can be computed at the end of each time step by: 

min(A;„Aj,) 
max{|^|,|f 1} 

271 where c is known as the Courant-Friedrichs-Lewy (CFL) number. To ensure the numerical stability of the time 

272 discretization scheme, we require that c < 1 . 
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273 4.3. Inversion subproblem 

274 Most of the demand on computing resources posed by QG models comes in the solution of the elliptic inversion 

275 subproblem ( [Miller 2QQ7| ). This is also true for our study. However, we take advantage of the simple square shape of 

276 our domain and utilize one of the fastest available techniques ( |Moin| [2001] [San and Staples] |2Q1 3 aj ), which is the FFT 

277 based direct inversion to solve the subproblem: 



Gi = RoVVi + ?(<A2-<Ai) 



22 = R0VV2 + 



Fr 



(47) 
(48) 



278 where Qi = qi-y and Q2 = - y- The impermeability boundary condition imposed as = suggests the use of 

279 a fast sine transform (an inverse transform) for each layer: 



A^,-1A^-1 



NxNy j 



:1 J-1 



nki 



nlj 



(49) 



N,-\ N,-\ 



(50) 



281 where Nx and Ny are the total number of grid points in x and y directions. Here the symbol hat is used to represent 

282 the corresponding Fourier coefficient of the physical grid data with a subscript pair /, 7, where / = 0, 1, ...Nx and 

283 J = 0, 1, ...A^3;. As a second step, we directly solve the subproblem in Fourier space: 



(^k,iQh,i - T^Qhj 



FrA 



( Fr Fr\ ' 
Fr, 



,l{(^k,l- ^ - ^) 



285 where 



o:k,l 



Ro 

A? 



2 cos 



Ro 



2 cos 



286 Finally, the streamfunction arrays for each layer are found by performing a forward sine transform: 



(51) 



(52) 



(53) 



^v-\^^v-n^ . (nkiX . l7ilj\ 



(54) 



ZY^ ./ TIKI \ , I TIL J 



Tiki 



Tllj 



288 The computational cost of this elliptic solver is 0(^NxNy log(A^;c) log(A/3;)). The FFT algorithm given by 

289 ( |1992| ) is used for forward and inverse sine transforms. 



(55) 



Press et al. 



290 5. Results 

291 The main goal of this section is to test the new AD model ([29|)-([3T|) in the numerical simulation of the two-layer 

292 QG model. We also investigate the sensitivity of the AD model with respect to the model parameters. It turns out 

293 that the most important modeling choice is the spatial filter employed in the AD procedure. We consider two spatial 

294 filters in conjunction with the AD model: the tridiagonal filter (Section [3?T]) and the diff'erential filter (Section [T2). 
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Table 2: Physical parameter sets used in the numerical experiments. 



Variable (unit) Experiment 1 Experiment 2 



L {KM) 


jUUU 


zuuu 


A 512^ A 512^ /J \ 

Ax = Ay^ (km) 


9.765625 


3.90625 


Hi (km) 


0.6 


1.0 


H2 (km) 


3.4 


4.0 


/o (s ') 


9.35 X 10 


9.35 X 10 


(m ^) 


1.75 X 10 


1.75 X 10 


pi (kgm ') 


1030 


1030 


/ / -2\ 

g (ms ^) 


0.02 


0.02 


To (Nm~^) 


0.1 


0.1 


y(s ') 


4 X 10"^ 


5 X 10"^ 


y (m s ) 


100 


50 


Sm (km) 


17.88 


14.19 


Ss (km) 


22.86 


2.86 


Sj (km) 


25.77 


31.56 


Rd (km) 


31.16 


42.79 


V(ms-^) 


0.0116 


0.0174 


L/V (year) 


13.64 


3.64 


LIRd 


160.5 


46.74 


Ro 


2.66 X 10-5 


2.49 X 10 


Fr 


0.073 


0.087 


cr 


4.57 X 10"^ 


1.43 X 10 


A 


4.57 X 10-^ 


3.57 X 10 


S 


0.15 


0.2 


Re 


580.97 


697.16 



We denote the resulting models AD-TF and AD-DF, respectively. To test the AD-TF and AD-DF models, we utilize 
two different parameter sets, corresponding to two physical oceanic settings: (i) Experiment 1 represents a large ocean 
basin with the physical parameters used by |Tanaka and Akitomo|(|2010|), (ii) E xperiment 2 represents a moderate ocean 
basin with the physical parameters used by 'Qz gokmen and Chassignet| ( |l99 8 ). In terms of the classification given by 



Berlolf and Mc Williams (1999), both sets of experiments lie under the chaotic regime. The physical parameters 



and corresponding dimensionless parameters are summarized in Table [2] All computations were carried out using a 



gfortran compiler on a Linux cluster system. The rest of the section is organized as follows. In Section 5.1 we present 



results from the direct numerical simulation (DNS) for the two settings. Experiment 1 and Experiment 2. Section 5.2 



presents results with the AD-TF model. Finally, Section 5.3 presents results with the AD-DF model. 

5.7. Direct numerical simulation 

We start by performing a DNS on a fine mesh of 512^ spatial resolution. We emphasize that the term DNS 
in this study is not meant to indicate that a fully detailed solution is being computed on the molecular viscosity 
scale, but instead refers to resolving the simulation down to the Munk scale via the specified lateral eddy viscosity 
parameterization. We also emphasize that the DNS results are given by the numerical solution of Eqs. ([29]) and 
( [3Q| ) with = S2 = 0. A statistically steady state solution is obtained after an initial transient spin-up process. 
Instantaneous contour plots for the potential vorticities in the upper and lower layers are shown in Fig. [3] and Fig. [4] 
for Experiment 1 and Experiment 2, respectively. The length scales in these two experiments are quite diff'erent. For 
example, the ratio of the basin length scale L to the Rossby deformation radius R^ is L/R^ = 160.5 for Experiment 1 
and L/Rd = 46.74 for Experiment 2. Therefore, the structure of the eastward jet formation on the western boundary 
for Experiment 1 is diff'erent from that of Experiment 2. This diff'erence becomes more obvious in the mean flow 
field. The results for time-averaged mean field data obtained from 2000 snapshots in the statistically steady state 
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Figure 3: Experiment 1: DNS results for (a) instantaneous potential vorticity for the upper layer, and (b) instantaneous potential vorticity for the 
lower layer. 



are given in Fig. [5] and Fig.|6] The results show strong western boundary currents with cyclonic (counter-clockwise 
rotating) subpolar gyres and anticy clonic (clockwise rotating) subtropical gyres producing a strong eastward jet in 
both experiments. However, the produced eastward jet formation in Experiment 2 shows swirling structure and almost 
reaches the eastern boundary of the basin. Compared to Experiment 1, the bottom lay er is more active in Experiment 2 . 
Since in Experiment 2 we used the same parame ters and b oundary conditions as in Qzgokmen and Chassignet ( 1998 ), 
the plot in Fig. [6] is similar to Fig. 2 in Ozgokmen and Chas signet] ( | 1998 ). Although in Experiment 1 we have used 
the same parameters as those used in Tanaka an d Akitomo| ( |20 1Q ), the boundary conditions we used are different from 
their boundary conditions: we used the slip boundary conditions, whereas they used the no-slip boundary conditions. 
Thus, the plot in Fig.|5]is different from the corresponding one in Tanaka and Akitomo (20T0|. 

To quantify the effect of the numerical discretization on the numerical results, we vary the grid resolution (NxXNy), 
the time step (AO, and the eddy viscosity coefficient (v) in the QG2 model. The following quantities are monitored. 
The first quantity is the time-averaged norm of the error of the potential vorticity, denoted as ||^/||, where the 
subscript / represents the layer index. The reference solution used in the computation of the error is the numerical 
approximation obtained at a grid resolution of 512^. The second quantity is the time-averaged basin-integrated kinetic 
energy, Ei, which is defined as 



Ei 



VP' 



Ei(t)dt, 



(56) 



Ei(t) : 



dxdy, 



(57) 



where, again, the subscript / represents the layer index and Ti = 6 and T2 = ^ are the temporal bounds for the 
averaging window. The integrand Ei(t) in ([56]) is the instantaneous basin integrated kinetic energy in each layer and is 
defined as 

First, we investigate the effect of the grid resolution on the numerical results. To this end, we fix the time ste] 
At = 2x 10"^ and vary the grid resolution, NxX Ny, and the eddy viscosity coefficient in the QG2 model, v. Table 
presents the time-averaged basin-integrated kinetic energy of the upper layer, Ei, defined in ( [56j ). This table shows 
that, for most grid resolutions, accurate results are obtained for the high values of the eddy viscosity coefficient, v. 
For the lowest values of y, however, the results are inaccurate at the lower grid resolutions, and relatively accurate at 
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Figure 4: Experiment 2: DNS results for (a) instantaneous potential vorticity for the upper layer, and (b) instantaneous potential vorticity for the 
lower layer. 



339 the higher grid resolutions. This behavior is natural, since when the grid spacing is larger than the Munk scale, the 

340 smallest scale are not resolved, thereby producing grid-scale variability in the solution, which degrades the accuracy. 

341 Table [4] presents the time-averaged norm of the error of the potential vorticity in the two layers, ||^i|| and ||^2||. 

342 This table shows that, as expected, the error decreases as the grid resolution increases. We note that this decrease in 

343 the error is faster for the high values of v. This behavior is similar to that observed in Table |3] Finally, the results in 

344 Table |4] and are also plotted in Fig.|7] This figure clearly shows that a second-order spatial accuracy is obtained for 

345 the high values of v, and a first-order spatial accuracy is obtained for the lowest values of v. However, one important 

346 thing to note in Fig.|7]is that the convergence approaches second-order when A^;^ = Ny aproaches 256. This is because 

347 the mininum number of grid points at which we can start to expect convergence is when the Munk scale is resolved, 

348 i.e., when Nx = Ny = 280 for Experiment 1. Therefore, we emphasize that the use of 512^ resolution should suffice 

349 for a DNS, although just barely. 

350 Next, we investigate the effect of the time step on the numerical results. To this end, we fix the the eddy viscosity 

351 coefficient, v = 100 m^s~^, and vary the grid resolution, Nx x Ny, and the time step. At. Table [s] presents the time- 

352 averaged basin-integrated kinetic energy of the two layers, Ei and E2. This table shows that, for a fixed spatial 

353 resolution, varying the time step does not yield a significant change in the numerical results. To perform the time- 

354 accuracy analysis for the third-order Runge-Kutta time integration scheme, we fix the grid resolution at 5 12^ and plot 

355 in Fig. [s] the norm of the error at ^ = 0.0075 for different time step sizes. The data obtained with = 5 x 10"^ 

356 is used as reference solution to compute the error norms. We present results at an early integration time to prevent 

357 the numerical instability that could appear later on as a result of the violation of the CFL criterion. The log-log plot 

358 in Fig. [8] clearly shows that the time integration scheme achieves the expected third-order temporal accuracy for both 

359 the upper and lower layers in Experiment 1 . To investigate the effects of the adaptive time discretization described in 

360 Section [4~2l we performed the same numerical experiments as those in Table [5j this time, however, using the adaptive 

361 time-stepping scheme with a fixed CFL number c = 0.95. This approach yielded the same qualitative results as those 

362 in TableE 

363 The above numerical studies quantify the effects of the numerical discretization described in Section |4] The 

364 following general conclusions can be drawn. The spatial discretization is optimal (second-order) for high values of 

365 the eddy viscosity coefficient, and is suboptimal (first-order) for the low values that we use in this study. The time 

366 discretization error appears to be dominated by the spatial discretization error. Indeed, for a fixed grid resolution. 
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Figure 5: Experiment 1: DNS results for (a) mean streamfunction contours for the upper layer, (b) mean potential vorticity contours for the upper 
layer, (c) mean streamfunction contours for the lower layer, and (d) mean potential vorticity contours for the lower layer. 
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Figure 6: Experiment 2: DNS results for (a) mean streamfunction contours for the upper layer, (b) mean potential vorticity contours for the upper 
layer, (c) mean streamfunction contours for the lower layer, and (d) mean potential vorticity contours for the lower layer. 
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Table 3: Experiment 1: Time-averaged basin-integrated kinetic energy of the upper layer, E\, for varying grid resolutions, Nx x Ny, varying eddy 
viscosity coefficients, v, and fixed time step At = 2x 10~^. 



N.XNy 


y = 100 


y = 200 


y = 400 


y = 800 


y = 1600 


y = 3200 




195.028 


200.188 


151.178 


89.139 


57.016 


36.500 


64^ 


103.787 


77.749 


59.083 


43.305 


33.567 


27.878 


128^ 


77.617 


63.618 


51.364 


42.545 


34.003 


27.661 


256^ 


79.478 


65.560 


52.646 


42.208 


34.764 


27.851 


512^ 


81.609 


66.084 


52.787 


42.051 


35.096 


27.921 



Table 4: Experiment 1: Time-averaged norm of the error of the potential vorticity in the two layers, ||^i|| and 11^2 1 L for varying grid resolutions, 
Nx X Ny, varying eddy viscosity coefficients, y, and fixed time step At = 2x 10~^. The reference solution used in the computation of the error is 
the numerical approximation obtained at a grid resolution of 512^. 



y = 100 



y = 400 



y = 3200 



qi 



I q\ 



qi 



32^ 
642 
128^ 
256^ 



1.2446E-1 
6.5465E-2 
2.9121E-2 
1.2296E-2 



1.8075E-2 
7.7261E-3 
3.3675E-3 
1.5355E-3 



1.4552E-1 
4.3220E-2 
1.3513E-2 
4.5496E-3 



2.0959E-2 
5.2517E-3 
1.7138E-3 
5.4868E-4 



4.7177E-2 
1.6356E-2 
4.7441E-3 
1.0002E-3 



7.2268E-3 
2.5420E-3 
7.4199E-4 
1.5671E-4 



Table 5: Experiment 1: 


Time-averaged basin-integrated kinetic energy of the two layers, E\ and E^, for varyin 


g grid resolutions, NxXNy, and fixed 


eddy viscosity coefficient, y = 100 m^^ ^ . 














= 1 X 10-5 




At = 2x 10-5 




= 4 X 10-5 






Ex 


El 


El 


E2 


El 


E2 


322 


198.862 


IA24 


195.028 


1.086 


196.293 


1.095 


64^ 


104.332 


0.874 


103.787 


0.876 


104.143 


0.875 


128^ 


78.210 


1.195 


77.617 


1.961 


77.768 


1.952 


256^ 


79.194 


2.532 


79.478 


2.523 


79.416 


2.538 


512^ 


81.277 


2.592 


81.609 


2.594 


80.996 


2.601 



changing the time step had only neghgible effects on the numerical results. Although it is hard to decouple the 
numerical effects from the LES modeling effects, the above numerical studies will serve as a guide in the subsequent 
interpretation of the LES results. Furthermore, a more detailed presentation of error estimates for the spatial and 
temporal schemes utilized here can be found in a recent study on two-dimensional decaying turbulence conducted by 
|San and Stales] ( [20T2| ). 



5.2. Approximate deconvolution model with the tridiagonal filter (AD-TF) 

To test the new AD-TF model ([3Q|-([3T]), we employ the standard LES methodology: We first run a DNS on a 
fine mesh (of 512^ spatial resolution). We then run on a much coarser mesh (of 32^ spatial resolution) an under- 
resolved numerical simulation (denoted in what follows as QG2c). We emphasize that QG2c does not employ any 
subfilter-scale model. Finally, we employ the new AD-TF model on the same coarse mesh utilized in QG2c (of 32^ 
spatial resolution). The criterion used in assessing the success of the new AD-TF model is its ability to produce more 
accurate (i.e., closer to the DNS data) results than those for QG2c, without a significant increase in computational 
time. Following San et al. ( 2011| ), in the AD-TF model, we use the tridiagonal filtering procedure with N = 5 and 



a = 0.25. To compare the DNS, the QG2c, and the AD-TF model, we utilize data that is time-averaged between t = 6 
and ^ = 8 by using 2000 snapshots of the field. Note that that this averaging period corresponds to 27.28 years for 
Experiment 1. 



17 



10-^ - 



O- 10-' - 



10-' - 






- v=100mV 


------ 


- v = 400mV^ 


— »„_„_, 


- v = 3200mV 



10-^ - 



32 



64 96 128 160 192 224256 

(a) 



10-^ 




V = 1 00 m^s" 

V = 400 m^s" 



v = 3200 mV 



32 



128 160 192 224256 



TV ( = # ) 
(b) 



Figure 7: Experiment 1: Log-log plot of the time-averaged norm of the error of the potential vorticity in the two layers, \\qi\\ and 11^2 IL for 
varying eddy viscosity coefficients v. The reference solution used in the computation of the error is the numerical approximation obtained at a grid 
resolution of 512^. 
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Figure 8: Experiment 1: Log-log plot of the norm of the error of the streamfunction, | 



I and I 



and potential vorticity, ||^i|| and \\q2\\, 



in the two layers at t = 0.0075, using the eddy viscosity coefficient v = 100 rn^s ^ and a resolution of 512^. The reference solution used in the 
computation of the error is the numerical approximation obtained by using the time step At = 5 x 10~^. 

For Experiment 1, we plot the mean streamfunction and potential vorticity contours in Figs.|9]and[T0j respectively. 
The new AD-TF model yields results that are significantly better than those corresponding to the under-resolved QG2c 
run. Similarly, we plot the mean streamfunction and potential vorticity contours in Figs. [TT] and [12] for Experiment 2. 
We note that the proposed AD-TF model yields again improved results by smoothing out the numerical oscillations 
present in the under-resolved QG2c simulations. We also note that the computational cost of the new AD-TF model is 
significantly lower than that of the DNS, and is comparable to the computational cost of the QG2c. Indeed, the CPU 



18 




Figure 9: Experiment 1: Time-averaged streamfunction contours for the upper layer: (a) DNS results at a resolution of 512^; (b) QG2c (under- 
resolved numerical simulation without any subfilter-scale model) results at a resolution of 32^; and (c) AD-TF results at a resolution of 32^. The 
contour layouts are identical. Note that the AD-TF results are significantly better than the QG2c results. 




Figure 10: Experiment 1: Time-averaged potential vorticity contours for the upper layer: (a) DNS results at a resolution of 512^; (b) QG2c 
(under-resolved numerical simulation without any subfilter-scale model) results at a resolution of 32^; and (c) AD-TF results at a resolution of 32^. 
The contour layouts are identical. Note that the AD-TF results are significantly better than the QG2c results. 



389 time is 1 19.6 hrs. for the DNS, 141.8 sees, for QG2c, and 174.7 sees, for the AD-TF model. The numerical results for 

390 both experiments clearly suggest that the the AD-TF model can provide relatively accurate results for under-resolved 

391 geophysical flows at a low computational cost. 

392 Although the AD-TF model performs well given the coarse mesh utilized, a natural question is whether we can 

393 increase its accuracy by using a finer mesh. Thus, we investigate the behavior of the AD-TF model for various 

394 resolutions: 512^, 256^, 128^, 64^, and 32^. Since similar conclusions hold for both experiments, we only discuss 

395 the results for Experiment 1 . The time-averaged streamfunction and potential vorticity contour plots for the upper 



400 
401 
402 
403 



396 layers are shown in Figs. 13 and 14 respectively. The qualitative results displayed in these figures are reinforced by 

397 the quantitative results in Table [6j which presents the time-averaged norm of the error of the potential vorticity 

398 in the two layers, ||^i|| and ||^2lL for fixed truncation order, N = 5, varying grid resolutions, A^;^ X and varying 



399 free parameter a. These results are also compared graphically in Fig. [15] The CPU time is 296 hrs for the DNS 
results, 48.5 hrs for the 256^ resolution, 4.1 hrs for the 128^ resolution, 0.34 hrs for the 64^ resolution, and 2.9 mins 
for the 32^ resolution. The main conclusion that can be drawn from the plots in Figs. [l3] [l4j [15] Table [6] and the 
computational efficiency study is that at the lowest resolutions the AD-TF model achieves a very high speed-up factor 
and an acceptable order of accuracy with respect to the DNS results (significantly higher than the accuracy of QG2c, 
404 i.e., the under-resolved numerical simulation without any subfilter-scale model). We also conclude that the AD-TF 
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Figure 11: Experiment 2: Time-averaged streamfunction contours for the upper layer: (a) DNS results at a resolution of 512^; (b) QG2c (under- 
resolved numerical simulation without any subfilter-scale model) results at a resolution of 32^; and (c) AD-TF results at a resolution of 32^. The 
contour layouts are identical. Note that the AD-TF results are significantly better than the QG2c results. 






Figure 12: Experiment 2: Time-averaged potential vorticity contours for the upper layer: (a) DNS results at a resolution of 512^; (b) QG2c 
(under-resolved numerical simulation without any subfilter-scale model) results at a resolution of 32^; and (c) AD-TF results at a resolution of 32^. 
The contour layouts are identical. Note that the AD-TF results are significantly better than the QG2c results. 
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405 model is consistent with the original set of equations, since the AD-TF results converge to the DNS results when the 

406 mesh size approaches zero. 

407 Finally, we perform a sensitivity study of the free smoothing parameter a and the order N in the AD-TF model. 

408 For comparison purposes, we also include results for QG2c (the under-resolved numerical simulation without any 

409 subfilter-scale model). In order to quantify the results of the AD-TF model, we compute the error norms with respect 

410 to the DNS results with a resolution of 512^. In both DNS and QG2c computations, the subfilter-scale term is set to 

411 zero: Sl= 3^ = 0. 



412 We start by investigating the sensitivity of the AD-TF model with respect to the parameter a. Table [6] and Fig. 15 



413 show that the sensitivity of the results to the free parameter a decreases with increasing mesh refinement. Indeed, 

414 at the coarsest resolution (i.e., 32^), the values a = 0.15, a = 0.25, and a = 0.35 yield practically indistinguishable 

415 results. The value a = 0.45 yields the most inaccurate results. At the 64^ resolution, the value a = 0.15 yields the 

416 best results, whereas the value a = 0.45 yields again the most inaccurate ones. At the 128^ and 256^ resolutions, the 

417 results are similar for all the values of a. In conclusion, the value a = 0.25 appears to be optimal, since it yields the 

418 best results at the 32^ resolution in Table[6]and Fig. 15 We note, however, that the values a = 0.15 and a = 0.35 yield 



similar results. We also note that, for low values of a, the AD-TF model performs better than QG2c at all resolutions. 

420 For higher values of a, the AD-TF model performs better than QG2c at the lowest resolution, but its accuracy starts 

421 to degrade at higher resolutions. 

422 Next, we investigate the sensitivity of the AD-TF model with respect to the order A^. To this end, in Table |7j we 

423 fix the parameter a = 0.25 and the grid resolution 32^, and present the time-averaged norm of the error of the 
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Figure 14: Experiment 1: Time-averaged potential vorticity for the upper layer: (a) DNS results at a resolution of 512^; (b) AD-TF results at a 
resolution of 512^; (c) AD-TF results at a resolution of 256^; (d) AD-TF results at a resolution of 128^; (e) AD-TF results at a resolution of 64^; 
and (f) AD-TF results at a resolution of 32^. The contour layouts are identical. Note: (i) the accuracy of the AD-TF results at the 128^ resolution; 
and (ii) the consistency of the AD-TF results with respect to the mesh size. 

424 streamfunction, \\i//i\\ and \\i/^2\l and potential vorticity in the two layers, ||^i|| and ||^2lL for varying orders N in the 



433 



440 



425 AD-TF model. These results are also compared graphically in Fig. [16] We note that, at this coarse resolution, the 

426 AD-TF model performs better than QG2c for all values of N. Based on the results in Table|7]and Fig.[T6j we conclude 

427 that the truncation order N = 3 is the optimal value for the AD-TF model. Indeed, increasing the value of N from 

428 1 to 3, results in a significant decrease in the error. For higher values of N, however, the decrease in the error is 

429 negligible. Since increasing the value of N implies more filtering operations in the computation of the subfilter- scale 

430 
431 



term and, thus, a higher computational time, the value N = 3 yields the best results in terms of combined accuracy 
and eflftciency. 



5.3. Approximate deconvolution model with the differential filter (AD-DF) 

Section[5]2]clearly showed that, for a fixed value of the EV coefficient v, the AD-TF model can provide an accurate 

434 approximation of the mean flow field on a mesh that is significantly coarser than that used in a DNS. Furthermore, it 

435 also showed that the AD-TF model is relatively insensitive with respect to changes in the smoothing parameter a used 

436 in the definition of the tridiagonal filter. A natural question is whether the AD model is sensitive with respect to other 

437 choices in the input parameters, such as the spatial filter. In this section, we numerically investigate the AD model 



438 equipped with a differential filter (given in Eq. ( [37] ) and discussed in Section ]T2] ) instead of the tridiagonal filter used 

439 in Section ]S!2] The resulting LES model is denoted as AD-DF. 



We start by performing a sensitivity study with respect to the model parameter A and the order N in the AD- 



DF model, similar to the analysis performed in Section 5.2 for the AD-TF model. For comparison purposes, we 
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Figure 15: Experiment 1: Log-log plot of the time-averaged norm of the error of the potential vorticity in the two layers, (a) ||^i|| and (b) \\q2\\, 
for varying parameter a in the AD-TF model. The results obtained with QG2c (the under-resolved numerical simulation without any subfilter-scale 
model) are also included for comparison purposes. The reference solution used in the computation of the error is the numerical approximation 
obtained at a grid resolution of 512^. 
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Figure 16: Experiment 1: The time-averaged norm of the error of the potential vorticities in the two layers, ||^i|| and ||^2lL for varying orders 
N in the AD-TF model, at a grid resolution of 32^, and for a fixed parameter a = 0.25. The results obtained with QG2c (the under-resolved 
numerical simulation without any subfilter-scale model) are also included (dashed lines) for comparison purposes. The reference solution used in 
the computation of the error is the numerical approximation obtained at a grid resolution of 512^. 
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Figure 17: Experiment 1: Log-log plot of the time-averaged norm of the error of the potential vorticity in the two layers, (a) ||^i|| and (b) 
11^211 for varying Helmholtz length A in the AD-DF model. The results obtained with QG2c (the under-resolved numerical simulation without any 
subfilter- scale model) are also included for comparison purposes. The reference solution used in the computation of the error is the numerical 
approximation obtained at a grid resolution of 512^. 




Figure 18: Experiment 1: The time-averaged norm of the error of the potential vorticity in the two layers, ||^i|| and ||^2lL for varying orders 
N in the AD-DF model, at a grid resolution of 32^, and for a fixed parameter A = 0.6/z. The results obtained with QG2c (the under-resolved 
numerical simulation without any subfilter-scale model) are also included (dashed lines) for comparison purposes. The reference solution used in 
the computation of the error is the numerical approximation obtained at a grid resolution of 512^. 



24 



Table 6: Experiment 1: Time- averaged norm of the error of the potential vorticity in the two layers, \\q\\\ and ||^2lL for varying grid resolutions, 
Nx X Ny, and varying the parameter a of the AD-TF model. The results obtained with QG2c (the under-resolved numerical simulation without 
any subfilter- scale model) are also included for comparison purposes. The reference solution used in the computation of the error is the numerical 
approximation obtained at a grid resolution of 512^. 

a = 0.15 a = 0.25 a = 0.35 a = 0.45 QG2, (5* = 0) 



qi II II qi II II qi II II qi II II qi II II qi II II qi II II qi II II qi II II qi 



32^ 


6.07E-2 


1.02E-2 


6.00E-2 


9.68E-3 


6.24E-2 


9.45E-3 


7.73E-2 


1.14E-2 


1.24E-1 


1.81E-2 


64^ 


4.36E-2 


6.42E-3 


5.06E-2 


7.50E-3 


5.69E-2 


7.82E-3 


7.08E-2 


8.48E-3 


6.55E-2 


7.73E-3 


128^ 


2.72E-2 


3.34E-3 


2.90E-2 


3.48E-3 


3.19E-2 


3.73E-3 


3.85E-2 


4.07E-3 


2.91E-2 


3.37E-3 


256^ 


1.07E-2 


1.39E-3 


l.lOE-2 


1.42E-3 


1.19E-2 


1.48E-3 


1.32E-2 


1.59E-3 


1.23E-2 


1.54E-3 



Table 7: Experiment 1: Time-averaged norm of the error of the streamfunction, ||i^i|| and ||<^2lL and potential vorticity in the two layers, ||^i|| 
and 11^211, for a fixed parameter or = 0.25, and varying orders N in the AD-TF model. The reference solution used in the computation of the error is 
the numerical approximation obtained at a grid resolution of 512^. 



Method (Njc x A^3;) 


II ^1 II 


II ^2 II 


II qi II 


II qi II 


QG2, (32^) 


2.2090E-1 


2.6845E-2 


1.2446E-1 


1.8075E-2 


AD-TF; = 1 (32^) 


2.0484E-1 


1.9669E-2 


1.1848E-1 


1.6965E-2 


AD-TF; = 2 (32^) 


1.4563E-1 


3.0302E-2 


8.0293E-2 


1.3159E-2 


AD-TF; = 3 (32^) 


1.1990E-1 


2.5521E-2 


6.1966E-2 


1.0139E-2 


AD-TF; = 4 (32^) 


1.1642E-1 


2.443 lE-2 


6.0171E-2 


9.7303E-3 


AD-TF; A^ = 5 (32^) 


1.1701E-1 


2.3525E-2 


5.9979E-2 


9.6823E-3 



Table 8: Experiment 1: Time-averaged norm of the error of the potential vorticity in the two layers, ||^i|| and ||^2lL for varying grid resolutions, 
Nx X Ny, and varying the parameter A of the AD-DF model. The results obtained with QG2c (the under-resolved numerical simulation without 
any subfilter-scale model) are also included for comparison purposes. The reference solution used in the computation of the error is the numerical 
approximation obtained at a grid resolution of 512^. 



N^XNy 


A = 0.4/1 




A = OM 




A = OM 




A = l.Oh 




QG2, (5* 
II qi II 


= 0) 


II qx II 


II qi II 


II qi II 


II qi II 


II qi II 


II qi II 


II qi II 


II qi II 


II qi II 


32^ 
642 
128^ 
256^ 


8.72E-2 
5.08E-2 
2.74E-2 
1.24E-2 


1.29E-2 
7.00E-3 
3.48E-3 
1.49E-3 


7.03E-2 
5.01E-2 
3.07E-2 
1.38E-2 


1.07E-2 
7.36E-3 
3.76E-3 
1.66E-3 


7.09E-2 
5.33E-2 
3.53E-2 
1.75E-2 


1.23E-2 
8.00E-3 
4.38E-3 
1.89E-3 


7.28E-2 
5.68E-2 
3.93E-2 
2.27E-2 


1.16E-2 
8.62E-3 
4.93E-3 
2.23E-3 


1.24E-1 
6.55E-2 
2.91E-2 
1.23E-2 


1.81E-2 
7.73E-3 
3.37E-3 
1.54E-3 



Table 9: Experiment 1: Time-averaged norm of the error of the streamfunction, and \\il/2\\, and potential vorticity in the two layers, ||^i|| 
and 11^211, for a fixed parameter A = 0.6h, and varying orders N in the AD-DF model. The reference solution used in the computation of the error is 
the numerical approximation obtained at a grid resolution of 512^. 



Method (A^;, x A^^;) 


II ^Ai II 


II 1/^2 II 


II qi II 


II qi II 


QG2, (32^) 


2.2090E-1 


2.6845E-2 


1.2446E-1 


1.8075E-2 


AD-DF; A^ = 1 (32^) 


3.3117E-1 


2.0303E-2 


1.9157E-1 


2.7892E-2 


AD-DF; A^ = 2 (32^) 


1.7354E-1 


2.0192E-2 


9.7860E-2 


1.4304E-2 


AD-DF; A^ = 3 (32^) 


1.4690E-1 


2.0361E-2 


8.0382E-2 


1.1957E-2 


AD-DF; A^ = 4 (32^) 


1.3297E-1 


2.0021E-2 


7.0474E-2 


1.0715E-2 


AD-DF; A^ = 5 (32^) 


1.3343E-1 


1.9745E-2 


7.0305E-2 


1.0723E-2 



also include results for QG2c (theunder-resolved numerical simulation without any subfilter-scale model). In order 
to quantify the results of the AD-DF model, we compute the error norms with respect to the DNS results having a 
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resolution of 512^. In both DNS and QG2c computations, the subfilter-scale term is set to zero: 5"* = 5"* = 0. 

We first investigate the sensitivity of the AD-DF model with respect to the Helmholtz length, A. To this end, in 
TablejsJ we fix the truncation order N = 5, and present the time-averaged norm of the error of the potential vorticity 
in the two layers, ||^i|| and ||^2ll, for varying grid resolutions, A^;^ X Ny, and varying Helmholtz length, A. These results 
are also compared graphically in Fig. 17 Table [s] and Fig. [It] yield the following conclusions. At the 32^ and 64^ 



resolutions, all the A values yield similar results. At the 128 and 256 resolutions, however, the values A = OAh and 
A = 0.6h yield the most accurate results; the values A = 0.8/z and A = I. Oh yield inaccurate results. In conclusion 
the value A = 0.6h appears to be optimal, since it yields the best results at the 32^ resolution in Table [s] and Fig. 
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We also note that, for the values A = OAh and A = 0.6h, the AD-DF model performs better than (or similar to) QG2c 
at all resolutions. For the values A = 0.8/z and A = I. Oh, however, the AD-DF model is more accurate than QG2c 
at low resolutions (32^ and 64^), but less accurate than QG2c at high resolutions (128^ and 256^). This behavior is 
natural, since, as explained in Section [^2] the higher values of A correspond to a higher level of numerical dissipation 
introduced by the DF. A higher level of dissipation is beneficial to the numerical simulations at low resolutions, since 
it models some of the subgrid-scale eff'ects. At higher resolutions, however, the subgrid-scale eff'ects become less 
important. In this case, the dissipation introduced by the DF should also decrease. This explains why, at higher 
resolutions, the lower values of A yield better results than the higher values of A. 

Next, we investigate the sensitivity of the AD-DF model with respect to the order N. To this end, in Table [9j we 
fix the filtering parameter at A = 0.6h and the grid resolution at 32^, and present the time-averaged norm of the 
error of the streamfunction, \\i/^i\\ and \\i/^2\l and potential vorticity in the two layers, ||^i|| and ||^2ll, for varying orders 
N in the AD-DF model. These results are also compared graphically in Fig. [18] Based on the results in Table [9] and 
Fig. 18 we conclude that the truncation order N = 4is the optimal value for the AD-DF model. Indeed, increasing the 
value of N from 1 to 4, results in a significant decrease in the error. For N = 5, however, the decrease in the error is 
negligible. Since increasing the value of N implies more filtering operations in the computation of the subfilter-scale 
term and, thus, a higher computational time, the value N = 4 yields the best results in terms of combined accuracy 
and eflficiency. 

The above sensitivity study clearly shows that, for a fixed value of the EV coefficient y, the AD-DF model can 
provide an accurate approximation of the mean flow field on a mesh that is significantly coarser than that used in 
a DNS. It was also shown that the diff'erential filter introduces a significant amount of numerical dissipation for 
higher values of A. The rest of the section is devoted to a careful numerical investigation of the amount of numerical 
dissipation in the AD-DF model by varying the EV coeflftcient v in the model. 

As mentioned in the introduction, the origin and modeling of the EV coefl&cient v in the QG models is a thorny 
issue (the "elephant in the room"). Indeed, Table [T] shows the wide range of values used for the EV coefficient v 
over the years. It is clear that no unique choice exists for y. Instead, the value used in numerical simulations is 
dictated by the available computational resources. To illustrate the importance of the particular value used for y in 
practical computations, we carried out several high-resolution 128^ numerical simulations for various EV coefficients 
y. Fig.[T9ja) shows the time series of the basin integrated kinetic energy for Experiment 1 for diff'erent values of y. The 
corresponding evolution of the maximum speed is also plotted in Fig.[T9jb), in which we convert the dimensionless 
velocity to its dimensional counterpart to get a better physical insight. As seen from Fig. 19 after an initial transient 
spin-up process, the system with y = 100 m^^"^ reaches a statistically steady state at an average maximum speed 
of 1.78 ms~^ (having an upper bound of 2.05 m^-"^ and a low er bound of 1.55 ms~^), which is close to the observed 
maximum zonal velocities of 2 ms~^ at 68°Ty (Dijkstra 2005 ). Thus, Fig. 19 illustrates the procedure used in choosing 
the EV coefficient y in practical computations with the QG model: The available computational resources dictate the 
numerical resolution that can be used; this, in turn, determines the EV coeflficient y that yields physical values for the 
computed flow fields (i.e., values that match those from observational data). Using higher or lower values for y can 



result in unphysical flow field data, as illustrated in Fig. 19 



In order to measure the amount of numerical dissipation in the AD-DF model with A = 2/z, we run this model 
with EV coefficients that span three orders of magnitude. The results for Re = 580.97 (y = 100 m^s~^), Re = 5809.7 
(y = 10 m^s~^), and Re = 58097 (y = 1 n?s~^) obtained with the AD-DF model are presented in Fig. 20 which shows 
the time histories of the basin integrated kinetic energy given by Eq. ([57]) for the AD-DF model (for all three Reynolds 
numbers). Results for the DNS and for the No-SFS run (the under-resolved numerical simulation without any LES 
model) for Re = 580.97 (y = 100 m^^"^) (the Reynolds number used in Section 5^) are also included for comparison 
purposes. As expected, No-SFS does yield a non-physical flow field with an unrealistically increasing energy level. 
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(a) (b) 

Figure 19: Experiment 1 : Numerical simulation at a 128^ resolution for different values of the EV coefficient v. (a) Time history of basin integrated 
kinetic energy given by Eq. {57} for the upper layer, and (b) time series of the maximum speed Vm in the field. Note the sensitivity of the results 
with respect to v. 



496 The kinetic energy of the AD-DF model for Re = 580.97 (v = 100 m^s~^), on the other hand, is significantly lower 

497 than the kinetic energy of the DNS. Thus, we conclude that the diff'erential filter in the AD-DF model yields too much 

498 numerical dissipation. Lowering the value of the eddy viscosity coefficient v alleviates this problem. Indeed, the AD- 

499 DF model with Re = 5809.7 (v = 10 n?s~^) produces the same level of kinetic energy as the DNS for the Re = 580.97 

500 (y = 100 m^^"^). Lowering even further the value of v results in a kinetic energy level that is unrealistically high. 

501 Based on the results in Fig. 20, we conclude that the diff'erential filter in the AD-DF model introduces a wider range 

502 of numerical dissipation in the model. 
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503 We have already seen in Fig. 20 that the AD-DF model can successfully run on coarse meshes at a lower eddy 

504 viscosity coefficient v. The transfer functions in Fig. [T] and Fig. [2] clearly suggest that the AD-DF model should 

505 provide more dissipation than the AD-TF model. The next natural step is to quantify how much eff'ective dissipation 

506 is provided by each model. We address this issue by performing numerical experiments for diff'erent values of the 

507 parameters a in the AD-TF model and X in the AD-DF model, and for various values of the eddy viscosity coefficient 

508 m^^"^ < y < 100 m^^"^ We note that, as expected, when the dissipation in the system is turned off' (i.e., v = m^^"^), 

509 both the DNS and QG2c computations do not reach a quasi- stationary energy level; this is indicated in Table [TOjby a 

510 dash symbol. The AD truncation order is fixed to A/^ = 5. The domain-integrated kinetic energy for the upper layer is 



511 presented in Table [T0| for diff'erent values of v, a, and A. The long time integrations are performed by using a coarse 

512 resolution of 32^ for all the runs. The results obtained by the QG2c, the under-resolved numerical simulation without 

513 any subfilter-scale model (i.e., 5* = 0) at the same resolution of 32^ and the DNS results obtained at a resolution 



of 512^ are also included for comparison purposes. Table 10 shows that the diff'erence in mean kinetic energy level 
between the DNS and QG2c is quantitatively high for all values of v. The reason is that the coarser resolution in 
QG2c does not eff'ectively resolve the Munk scale. Using the same coarse resolution of 32^, the AD-TF model with 
0.05 < a < 0.45 predicts a more accurate energy level. Decreasing a adds more numerical dissipation, which results 
in a decrease of the predicted energy level. For the higher values of a, the accuracy of the AD-TF model considerably 
degrades. Thus, when the parameter a is between a = 0.25 and a = 0.35, the AD-TF model yields the most accurate 

520 results. At the same coarse resolution of 32^, the AD-DF model with 0.4/z < A< 2.0h predicts a more accurate energy 

521 level than QG2c. Decreasing A adds a low level of numerical dissipation. Increasing A provides a significant amount of 

522 numerical dissipation and decreases the predicted energy level. The most accurate results are obtained by the AD-DF 
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Figure 20: Experiment 1: Time histories of basin integrated kinetic energy for the upper layer obtained by the AD-DF model for Re = 580.97 
(y = 100 m^^-i), Re = 5809.7 (y = 10 m^s-^), and Re = 58097 (y = 1 m^^-i). Results for the DNS and the No-SFS (under-resolved numerical 
simulation without any LES model) are also included. Note that the AD-DF introduces numerical dissipation and produces the same mean flow 
field as the DNS on a coarser mesh and for a lower eddy viscosity coefficient y. 
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model when the free parameter A Hes between A = OAh and A = 0.6h for all values of v. Overall, we conclude that, for 
a fixed value of the eddy viscosity coeflftcient y, we can obtain results close to the DNS results by tuning the modeling 



parameters a and J3 in the AD-TF and AD-DF models, respectively. Furthermore, the results in Table 10 show that 
the kinetic energy level predicted by the DNS for a given value of v can be predicted by the AD-TF and AD-DF 
models with a lower value of v when the model parameters a and j3 are appropriately chosen. Thus, as expected (see 
the transfer functions in Fig. [T] and Fig. [2]), the AD-TF and AD-DF models do provide numerical dissipation to the 
system. We note, however, that using the DF and TF without using the AD procedure does not provide a significant 
amount of numerical dissipation. Indeed, running an under-resolved numerical simulation with v = Om^s~^ and using 
the DF to smooth out the potential vorticity and streamfunction values after each time-step yielded inaccurate results. 
Finally, we emphasize that the numerical results in Table [T0| should not interpreted as an argument for the superiority 
of the AD-LES models over standard EV models. Instead, they simply show that the same kinetic energy level can 
be predicted in two diff'erent ways: by adjusting the eddy viscosity coefficient or by adjusting the parameters a and jS 



in the AD-TF and AD-DF models. For completeness, in Table 1 1 we repeat the same numerical experiments as those 
displayed in Table 10 but for a moderate resolution of 128^. The same qualitative conclusions as those above can be 



drawn, except that the diff'erence between DNS and QG2c in Table 1 1 is smaller due to the fairly well resolved Munk 
scales at this resolution. 



539 6. Conclusions 

540 A new approximate deconvolution large eddy simulation (AD-LES) model for the two-layer quasigeostrophic 

541 equations, a standard prototype of more realistic wind-driven ocean circulation, was introduced. Two diff'erent ocean 

542 settings with eastward jet formations of diff'erent strengths were considered. Two variants of the AD-LES model 

543 were proposed: one with a tridiagonal filter (AD-TF), and the other with a diff'erential filter (AD-DF). Both the 

544 AD-TF and the AD-DF models yielded accurate solutions, with physically relevant energy levels and realistic mean 

545 streamfunction and potential vorticity contour plots. A quantitative analysis of the eff'ects of using AD-TF and AD- 

546 DF on the QG2 model was presented. The two models also dramatically decreased the computational cost of the 

547 corresponding high-resolution numerical simulation, by using a mesh significantly coarser than the Munk scale. We 

548 emphasize that the AD procedure plays an essential role in the success of the AD-LES modeling strategy. Indeed, the 
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Table 10: Experiment 1: Time-averaged basin-integrated kinetic energy of the upper layer, E\, for varying modeling parameters, a and A, varying 
eddy viscosity coefficients, v, and a fixed resolution of 32^. The reference DNS solution obtained at a grid resolution of 512^ and the solution 
obtained by the QG2c without using any subfilter-scale model at a resolution of 32^ are also included for comparison purposes. The dash represents 
a nonstationary field. 



Method 




y = nt^s ^ 


y = 1 m^s ^ 


y = 5 ni^s ^ 


y = 10 m^s~^ 


y = 50 m^s~^ 


y = 100 m^s~^ 


AD-TF: a 


= 0.05 


68.744 


61.655 


52.597 


48.790 


33.884 


27.695 


AD-TF: a 


= 0.15 


69.620 


67.875 


60.054 


55.825 


41.881 


38.532 


AD-TF: a 


= 0.25 


85.855 


83.992 


78.974 


71.994 


53.641 


48.478 


AD-TF: a 


= 0.35 


107.444 


106.521 


94.917 


87.655 


83.846 


74.202 


AD-TF: a 


= 0.45 


369.895 


323.502 


245.234 


210.642 


145.469 


119.488 


AD-DF: A 


= OAh 


220.331 


212.155 


164.394 


139.083 


96.006 


84.913 


AD-DF: A 


= 0.6h 


135.768 


117.392 


90.414 


77.151 


50.566 


42.623 


AD-DF: A 


= 0M 


85.939 


86.124 


64.816 


51.566 


31.910 


26.804 


AD-DF: A 


= I. Oh 


71.628 


75.752 


49.740 


40.499 


23.448 


18.161 


AD-DF: A 


= 2.0/1 


49.906 


43.736 


27.501 


21.936 


9.556 


7.164 


QG2, (s; 


= 0) 




397.121 


333.205 


279.101 


210.718 


195.028 


DNS (5* = 


= 0) 




128.698 


121.196 


117.014 


96.593 


81.609 



Table 11: Experiment 1: Time-averaged basin-integrated kinetic energy of the upper layer, E\, for varying modeling parameters, a and A, varying 
eddy viscosity coefficients, v, and a fixed resolution of 128^. The reference DNS solution obtained at a grid resolution of 512^ and the solution 
obtained by the QG2c without using any subfilter-scale model at a resolution of 128^ are also included for comparison purposes. The dash represents 
a nonstationary field. 



Method 




y = ftp's ^ 


y = 1 rrP's ^ 


y = 5 nP's ^ 


y = 10 m^^-^ 


y = 50 m^s-^ 


y = 100 m^s-^ 


AD-TF: a 


= 0.05 


101.317 


95.878 


84.092 


76.945 


60.018 


53.374 


AD-TF: a 


= 0.15 


129.586 


124.196 


105.293 


94.670 


69.004 


59.290 


AD-TF: a 


= 0.25 


164.488 


159.169 


135.967 


118.318 


79.629 


66.723 


AD-TF: a 


= 0.35 


217.861 


218.466 


173.868 


147.645 


90.977 


74.544 


AD-TF: a 


= 0.45 


336.068 


393.889 


238.706 


190.372 


105.839 


83.325 


AD-DF: A 


= OAh 


240.249 


187.409 


125.260 


102.714 


69.399 


60.189 


AD-DF: A 


= 0.6h 


128.823 


105.917 


78.761 


68.464 


54.023 


48.134 


AD-DF: A 


= 0M 


80.713 


67.224 


52.469 


47.975 


41.991 


37.782 


AD-DF: A 


= I. Oh 


57.065 


47.336 


38.009 


35.994 


32.924 


30.316 


AD-DF: A 


= 2.0h 


26.376 


19.411 


15.607 


15.116 


14.149 


14.166 


QG2, (SI 


= 0) 




442.829 


244.128 


179.348 


96.646 


77.617 


DNS (S* = 


= 0) 




128.698 


121.196 


117.014 


96.593 


81.609 



549 underresolved numerical simulations without AD modeling on the same coarse mesh as that employed by the AD-LES 

550 models produced inaccurate results. The AD-TF and AD-DF models, however, had different behaviors in terms of the 

551 numerical dissipation added to the system. In fact, our numerical results showed that the AD-TF and AD-DF models 

552 can be employed successfully on meshes that are significantly coarser than the Munk scale and with an eddy viscosity 

553 coefficient that is dramatically lower than that used in the original two-layer quasigeostrophic equations by tuning the 

554 free parameters a and A appropriately. We emphasize that the tuning of the AD-LES model parameters is essential 

555 in obtaining accurate results. We also note that this paper does not claim the superiority of the AD-LES method 

556 over other eddy viscosity type closure approaches since the underlying quasigeostrophic equations utilize an intrinsic 

557 eddy viscosity coefficient to account for large scale dissipation. With this in mind, we also highlight that assessments 

558 and evaluations of various turbulence closure models for large eddy simulations of realistic oceanic basins are highly 

559 desirable, a topic we intend to further investigate in a future study. 
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